Linear Time-Invariant (LTI) Systems
Table of Contents
A system $\mathcal{H}$ is linear if it satisfies the following two properties:
1) Scaling $$\mathcal{H}\{ \alpha x \} = \alpha \mathcal{H} \quad \forall \, \alpha \in \mathbb{C}$$

2) Additivity
$$ \text{If}\,\, y_1 = \mathcal{H} \{ x_1 \} \,\, \text{and}\,\, y_2 = \mathcal{H} \{ x_2 \} \, \, \text{then } \mathcal{H} \{x_1 + x_2\} = y_1 + y_2$$
A system $\mathcal{H}$ processing infinite-length signals is time-invariant (shift-invariant) if a time shift of the input signal creates a corresponding time shift in the output signal

We usally consider Linear Time-Invariant (LTI) systems.
Convolution is defined as the integral of the product of the two functions after one is reversed and shifted
$$y[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] = x[n]*h[n]$$Output $y[n]$ came out by convolution of input $x[n]$ and system $h[n]$
%%html
<center><iframe src="https://www.youtube.com/embed/Ma0YONjMZLI?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
For Infinite-Length Signals
$$y[n] = x[n]*h[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] $$
$$ = \cdots \,h[n+2]\,x[-2] + h[n+1]\,x[-1] + h[n]\,x[0] + h[n-1]\,x[1] + h[n-2]\,x[2]\, \cdots $$
It is inner product of $h$ vectors and $x$
Convolution is product of matrix $H$ and $x$

(impulse) : $\delta[n] = \begin{cases} 1 \quad n = 0\\ 0 \quad \text{otherwise}\\ \end{cases}$

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import IPython.display as ipd
import librosa.display
from scipy import signal
# Convolve a rectangular pulse with itself
# pulse
N = 8
n = np.arange(0, N)
x = [0,0,0,1,1,1,0,0]
# Convolve pulse with itself
y = np.convolve(x, x)
# plot
plt.figure(figsize=(10, 8))
plt.subplot(2,1,1)
plt.stem(x, markersize=4)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse')
plt.subplot(2,1,2)
plt.stem(y, markersize=4)
plt.xlim([0, 14])
plt.ylim([-0.2, 3.2])
plt.title('pulse * pulse')
plt.show()

moving average (MA) filter
low-pass filter in time domain
# piecewise smooth signal
N = 50
n = np.arange(0, N)
V = np.hstack([np.ones([1, np.int(N/2)]), -np.ones([1, np.int(N/2)])])
x = np.sin(np.pi/N*n)*V
xn = x + 0.1*np.random.randn(1, N)
# construct moving average filter impulse response of length M
M = 5
h = np.zeros([1, M])
h[0:M] = 1/M
# convolve noisy signal with impulse response
y = signal.convolve(xn, h)
# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')
plt.subplot(2,2,2)
plt.stem(xn.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')
plt.subplot(2,2,3)
plt.stem(h.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')
plt.subplot(2,2,4)
plt.stem(y.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()
# haar wavelet edge detector
M = 4
h = np.zeros([1, M])
h[0,0:2] = -1/M
h[0,2:4] = 1/M
# convolve noisy signal with impulse response
y = signal.convolve(xn, h)
# plot
plt.figure(figsize=(12,8))
plt.subplot(2,2,1)
plt.stem(x.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal')
plt.subplot(2,2,2)
plt.stem(xn.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('piecewise smooth signal + noise')
plt.subplot(2,2,3)
plt.stem(h.T, 'b', 'w')
plt.xlim([-1, N])
plt.ylim([-1.5, 1.5])
plt.title('impulse response')
plt.subplot(2,2,4)
plt.stem(y.T, 'b', 'w')
plt.xlim([0, N])
plt.ylim([-1.5, 1.5])
plt.title('convoluted output')
plt.show()
%%html
<center><iframe src="https://www.youtube.com/embed/I3isHB9Iy7E?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
%%html
<center><iframe src="https://www.youtube.com/embed/ZW6pw89cuEA?rel=0"
width="560" height="315" frameborder="0" allowfullscreen></iframe></center>
x, sr = librosa.load('./data_files/violin_origional.wav')
x = x/max(x) # normalized
ipd.Audio('./data_files/violin_origional.wav', rate=sr) # play a wave file with sampling rate sr
# plot wave file
plt.figure(figsize=(10, 4))
librosa.display.waveplot(x, sr=sr)
plt.xlabel('time in sec')
plt.show()
# impulse response in a closed room (by gunshot)
h, sr = librosa.load('./data_files/gunshot.wav')
ipd.Audio('./data_files/gunshot.wav', rate=sr) # play a wave file with sampling rate sr
# plot wave file
plt.figure(figsize=(10, 4))
librosa.display.waveplot(h, sr=sr)
plt.xlabel('time in sec')
plt.show()
y = signal.convolve(x, h)
y = y/max(y)
# plot
plt.figure(figsize=(10,8))
plt.subplot(2,1,1)
librosa.display.waveplot(x, sr=sr)
plt.xlim([0, 6])
plt.title('original')
plt.subplot(2,1,2)
librosa.display.waveplot(y, sr=sr)
plt.xlim([0, 6])
plt.title('convoluted')
plt.show()
ipd.Audio(y, rate=sr) # play a wave file with sampling rate sr
import cv2
import skimage
# noise image
img = cv2.imread('.\image_files\lena_sigma25.png', 0)
print(img.shape)
# plot
plt.figure(figsize = (6,6))
plt.imshow(img, cmap = 'gray')
plt.axis('off')
plt.show()
# smoothing or noise reduction
M = np.ones([3,3])/9
img_conv = signal.convolve2d(img, M, 'valid')
# plot
plt.figure(figsize = (12,6))
plt.subplot(1,2,1)
plt.imshow(img, cmap = 'gray')
plt.title('Noisy Image')
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(img_conv, cmap = 'gray')
plt.title('Smoothed Image')
plt.axis('off')
plt.show()
# original image
img = cv2.imread('.\image_files\lena.png', 0)
# plot
plt.figure(figsize = (6,6))
plt.imshow(img, cmap = 'gray')
plt.axis('off')
plt.show()
# guess what kind of image will be produced after convolution
Mv = [[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]]
Mv2 = [[1, 0, -1],
[2, 0, -2],
[1, 0, -1]]
img_conv = signal.convolve2d(img, Mv, 'valid')
img_conv2 = signal.convolve2d(img, Mv2, 'valid')
# plot
plt.figure(figsize = (12,6))
plt.subplot(1,2,1)
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.title('Mv')
plt.axis('off')
plt.subplot(1,2,2)
plt.imshow(img_conv2, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv2))
plt.title('Mv2')
plt.axis('off')
plt.show()
# guess what kind of image will be produced after convolution
Mh = [[-1, -2, -1],
[0, 0, 0],
[1, 2, 1]]
img_conv = signal.convolve2d(img, Mh, 'valid')
# plot
plt.figure(figsize = (6,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')
plt.show()
M = np.array(Mv + Mh)/2
img_conv = signal.convolve2d(img, M, 'valid')
# plot
plt.figure(figsize = (6,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')
plt.show()
M = [[0, -1, 0],
[-1, 4, -1],
[0, -1, 0]]
img_conv = signal.convolve2d(img, M, 'valid')
# plot
plt.figure(figsize = (6,6))
plt.imshow(img_conv, cmap = 'gray', vmin = 0, vmax = np.amax(img_conv))
plt.axis('off')
plt.show()
$$ y[n] = \text{median}\{x[n-k],\cdots,x[n+k]\}$$
There are nonlinear neighborhood operations that can be perfromed for the purpose of noise reduction that can do a better job of preserving edges than simple smoothing filters.
Median filters can do an excellent job of rejecting certain types of noise, in particular, "shot" or impulse noise (outlier in a time series) in which some individual pixels or signals have extreme values.
x, sr = librosa.load('./data_files/violin_origional.wav')
x = x/np.max(x) # normalized
ipd.Audio(x, rate=sr) # play a wave file with sampling rate sr
# generate an audio signal with a salt and pepper noise
shot_noise = np.zeros(x.shape)
shot_noise = skimage.util.random_noise(shot_noise, mode='s&p', seed=4)
x_noise = x + shot_noise - np.mean(shot_noise)
ipd.Audio(x_noise, rate=sr)
# apply a linear low-pass filter
h = np.array([1/3, 1/3, 1/3])
x_avg = signal.convolve(h, x_noise)
ipd.Audio(x_avg, rate=sr)
# does not work very well
# apply a nonlinear filter
x_median = signal.medfilt(x_noise, 7)
ipd.Audio(x_median, rate=sr)
# WOW !!!
# plot
plt.figure(figsize=(10,16))
plt.subplot(4,1,1)
plt.plot(np.arange(0, len(x))/sr, x)
plt.ylim([-1, 1])
plt.subplot(4,1,2)
plt.plot(np.arange(0, len(x))/sr, x_noise)
plt.ylim([-1, 1])
plt.subplot(4,1,3)
plt.plot(np.arange(0, len(x)+2)/sr, x_avg)
plt.ylim([-1, 1])
plt.subplot(4,1,4)
plt.plot(np.arange(0, len(x))/sr, x_median)
plt.ylim([-1, 1])
plt.show()
img = cv2.imread('.\image_files\lena.png', 0)
img_shot = skimage.util.random_noise(img, mode = 's&p')
# plot
plt.figure(figsize = (12,6))
plt.subplot(1,2,1) # original image
plt.imshow(img, cmap = 'gray', vmin = 0, vmax = np.amax(img))
plt.axis('off')
plt.subplot(1,2,2) # Degraded image with Salt & Pepper noise
plt.imshow(img_shot, cmap = 'gray', vmin = 0, vmax = np.amax(img_shot))
plt.axis('off')
plt.show()
# see the performance
M = np.ones([3,3])/9
img_avg = signal.convolve2d(M, img_shot, 'valid')
# plot
plt.figure(figsize = (6,6))
plt.imshow(img_avg, cmap = 'gray', vmin = 0, vmax = np.amax(img_avg))
plt.axis('off')
plt.show()
# Median Filter
img_median = signal.medfilt(img_shot,5)
# plot
plt.figure(figsize = (6,6))
plt.imshow(img_median, cmap = 'gray', vmin = 0, vmax = np.amax(img_median))
plt.axis('off')
plt.show()
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')